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Abstract 

We consider a mathematical model comprising of four coupled ordinary differential equations (ODEs) 
for studying the hepatitis C (HCV) viral dynamics. The model embodies the efficacies of a combina- 
tion therapy of interferon and ribavirin. A condition for the stability of the uninfected and the infected 
steady states is presented. A large number of sample points for the model parameters (which were phys- 
iologically feasible) were generated using Latin hypercube sampling. Analysis of our simulated values 
indicated approximately 24% cases as having an uninfected steady state. Statistical tests like the x^-test 
and the Spearman's test were also done on the sample values. The results of these tests indicate a dis- 
tinctly differently distribution of certain parameter values and not in case of others, vis-a-vis, the stability 
of the uninfected and the infected steady states. 
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1 Introduction 

Hepatitis C, which is an infectious disease caused by the hepatitis C virus (HCV), has a widespread 
global prevalence. It is estimated to have infected 170 million individuals worldwide mill- HCV is usually 
transmitted through blood contact with an infected person. In places like United States and Europe, the 
primary mode of transmission of HCV is through injected drug usage d. In India however, the lack of 
effective and reliable anti-HCV screening amongst blood donors is a major source of infection New 
HCV infections (through intravenous drug usage or otherwise) can be primarily classified as acute and 
chronic |[T1. While the acute cases in general do not have major detrimental implications and the virus 
is cleared, the ramifications of chronic cases are far reaching |[T1. It could lead to the occun^ence of liver 
cirrhosis and may eventually lead to hepatocellular carcinoma (HCC). About 50 — 80% of HCV infections 
are chronic in nature ||5|. Of these 10 — 20% develop liver cirrhosis of which about 5% are likely to be 
afflicted with HCC 0. 

The current treatment for HCV infection involves the combination therapy of pegylated interferon (IFN) 
and ribavirin Q, which yields long term response in only 50% of the cases. There is very little therapeuitic 
alternative for the cases of non-responders. Feld et al. IH dwell on the mechanism of action of combina- 
tion treatment of IFN and ribavirin in HCV infected patients. They report limited success (6 — 12% for 
a six-month treatment and 16 — 20% for a one-year treatment) in case of an IFN based monotreatment. 
The combination therapy however resulted in significant improvement (more than 50%) in response rates. 
Perelson et al. Q in their review article, discuss the role of this combination therapy for HCV infection, 
taking into account the fall in the efficacy levels of the drugs between dosing intervals. 
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Several mathematical models have been proposed for the study of hepatitis C viral dynamics. One of 
the earhest mathematical models was proposed by Neumann et al. ||8l, wherein the dynamics of HCV and 
the effect of interferon-a-2b were studied in vivo. The quantitative representation of the dynamics involved 
the incorporation of some aspects of the earlier successful models for HIV and HBV El HI. The model 
involved three coupled ODEs, where the key factors were identified as uninfected hepatocytes, productively 
infected hepatocytes and free HCV virions. The model assumed the growth of the uninfected hepatocytes 
at a constant rate accompanied by a natural death rate. In addition, the HCV was modeled as infecting the 
hepatocytes, which also had a natural death rate. The infected hepatocytes in turn abetted the growth of HCV. 
In addition, a term for the efficacy of IFN was incorporated to reduce the rate of infection of hepatocytes and 
to block the production of HCV from infected hepatocytes. The model exhibited a decline in the levels of 
infected hepatocytes and HCV accompanied by an increase in the level of uninfected hepatocytes. A single 
phasic decline was observed in case the efficacy of IFN in blocking the production of virions was taken to 
be zero, which is inconsistent with clinical observations EKH. Otherwise it showed a biphasic decline in 
viral load which is more realistic from the biomedical point of view. The model indicated that the key role 
of IFN is in blocking the production of virions from infected cells as compared to blocking infection of 
hepatocytes, which is minimal HHH. 

Dixit et al. Q extended the work of Neumann et al. HI in the sense of explicitly including the action 
of ribavirin. In this model El 13 the virion population is divided into infectious and non-infectious. Dixit et 
al. 191 assume that ribavirin (either on its own or in conjunction with IFN) renders a fraction of the newly 
produced HCV non-infectious. Another assumption of the model (in contrast to Neumann et al. HI) was 
that the role of IFN in blocking the production of infected hepatocytes is not significant. The model predicts 
that ribavirin does not have an impact on the IFN induced first phase decline |21. If the efficacy of IFN is 
large enough, then the second phase decline is not significantly affected by ribavirin. However, if the IFN 
efficacy is much smaller than 1, then the impact of ribavirin in the second phase decline is more profound. 
A plausible explanation for this is that a high IFN efficacy results in low HCV levels, which in turn reduces 
the role that ribavirin plays in rendering the virions as non-infectious. The model was able to explain why 
ribavirin enhances the second phase of decline in some cases and not in others El- 

Dahari et al. ifTOl further advanced the work of Neumann et al. HI and Dixit et al. El. Their model 
ifTOl incorporated the proliferation of both infected and uninfected hepatocytes which was not considered in 
the earlier models. They include density dependent proliferation for both infected and uninfected hepato- 
cytes, which restricted the growth of liver to a maximum possible size. This model was able to explain the 
limitations of the model of Dixit et al. lO, which could not explain the non-response of some patients and 
the triphasic decline patterns. The model defined a critical threshold efficacy below which there cannot be 
sustained long term viral load decline ||9l[T0l|. In then case when the efficacy is above the threshold one ob- 
serves a decline in the viral load, which could be biphasic or triphasic. The triphasic decay can be attributed 
to the homeostatic liver regeneration. 



The model under consideration is a fusion of the earlier models |l8l|9l[T0l[TTl| and is given by the following 
system of four coupled ODEs, 
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Here T{t) and I{t) are the number of uninfected and infected hepatocytes respectively while Vi{t) and 
VNi{t) are the number of infectious and non-infectious virions (corresponding to HCV RNA genome equiv- 
alents) respectively. The assumption in the model is that the hepatocytes are being produced at a constant 
rate s and have a natural death rate of d, while proliferating at a rate r with Tmax being the maximum possi- 
ble hepatocyte (both uninfected and infected) population level. The hepatocytes are assumed to be infected 
by the virions at a rate /3. The proliferation of infected hepatocytes is also assumed to take place at the 
rate r with a natural death rate 5. In the absence of any kind of treatment, the infected hepatocytes produce 
infectious virions at a rate p. The administration of IFN lowers this production by a factor of (1 — Cp), where 
ep is the efficacy of IFN. Finally, the model incorporates the efficacy p of ribavirin which renders a fraction 
of the infectious virions non-infectious, with a death rate c for both these populations. The model under 
consideration admits two steady states, viz- the uninfected and the infected steady states, as given below 
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Under the physiological conditions of r > d and s < dT^asji it can be shown ifTOll that the condition for 
stability of the uninfected steady state is, 
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p/3r(«)T„,ax 

while that of the infected steady state is, 
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Thus there is a transcritical bifurcation at, 
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3 Latin Hypercube Sampling 



In this article we use the method of Latin hypercube sampling in order to generate a collection of param- 
eter values from a multivariate normal distribution with a specified mean vector l-i and covariance matrix S. 
The theoretical and computational aspects of this method can be found in the pioneering work of McKay et 
al. |[T2l which was further advanced by Stein |[T3l . Latin hypercube sampling is a multidimensional version 
of unidimensional stratification. llT4l . Latin hypercube sampling is used to generate parameter values from a 
given distribution. In this method, random variates are generated from a d-dimensional uniform distribution 
over the hypercube [0, 1)'^. Random variates from other d- dimensional distributions can be generated from 
the uniform variates using the inverse transform method |[T4l . Stratification in one-dimensional uniform 
distribution can be achieved by dividing the interval into K strata. The extension of this idea to multidimen- 
sional stratification is not trivial as the sample size of K'^ is computationally expensive even for moderate 
values of d (since K has to be large enough). We present a brief outline of this method (after Glasserman 

HI). 

To begin with, we generate random variates V-^'^^ from a uniform distribution over [{j — 1)/K, j / K) 
for i = 1,2, ... ,d and a set of independent permutations ai, . . . ,ad over the set {1, 2, ... , K}{from K\ 

possibilities). The vectors {yi\ V^2^"'^ • • • ' ^P^) (for 3 — 1) 2, . . . , K) represent uniformly distributed 

points over the hypercube [0, l)'^. We can then set ^ V"'^^'' which gives us a stratified sample over the 
hypercube. An obvious way to achieve this is to set, 

y{j) ^ + \ i = l,...,d, j = l,...,K, and U^^^ ~ Uniform [0, 1) 

In generating this construction, one of the crucial assumptions is the independence of marginals, for in- 
troducing a correlation between parameter distributions might alter the stratification properties of Latin 
hypercube sampling. This is very much evident in generating a sample from a normal distribution which 
does not have a diagonal covariance matrix. The values so generated will not generally be stratified. In order 
to generate from (say) a normal distribution with mean vector and diagonal covariance matrix 

S = {^ij)ij=i we set, 

Z? = %!^., {v}'^),^ = l,■■■,d,j = l,...,K 
where ^^^.e^^ is the one-dimensional cumulative normal distribution with mean fii and variance Sjj. 



4 Results and Discussion 

We implemented the Latin hypercube sampling procedure describe above in MatLab^"^ and generated 
K = 2^^ sample points for all the parameter values. The mean vector /i was taken to be the values 
given in Table [T] i.e, fi = {s,d,p, l3,c,6,r,T^gj^)^ . The diagonal covariance matrix S was taken as 
diag (s^, ci^,p^, c^, J^, r^, T^j^j^). Once the K = 2^^ sets of parameter values were generated, they 
were tested for positivity as well as the physiological conditions (r > d and s < dTmax) as given in Section 
2. This test left us with 8213 sample parameter sets that were feasible. 

A number of tests were performed to understand how the sample parameter values correlate. Firstly, a 
test was performed on the 8213 accepted sample points against the frequency of positive samples from 
their respective intended distributions under the goodness of fit null hypothesis. The null hypothesis (set 
of sample points which are biomedically feasible do not deviate from the intended distribution), was not 
rejected for any of the sample parameters. On the other hand, a test done under the independence null 
hypothesis rejects this hypothesis for all parameters with respect to the significance level chosen (P < 0.05). 
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Since, physiological conditions add "restrictions" to the intended distribution, we can surmise from these 
tests that the accepted values do not deviate from this distribution. 

The accepted sample points were categorized into three parts depending on whether the value 

p/3r(")r^ax 

was greater than 1, less than or was in the interval [0, 1]. Notice that, when < Cr* < 1, both steady 
states can be stable on two disjoint sets of feasible Ep and p values (due to the transcritical bifurcation). The 
number of such sample points, which emerged from our simulations, was 6262 (76.25%). On the other 
hand, when Cr* > 1 or Cr* < 0, only one of the two steady states is stable. The sample parameter values, 
which gave the Cr* values to be greater than 1 can only lead to the stability of the uninfected steady state, 
since in this case the stability condition for the uninfected steady state is satisfied for all possible values of 
ep or p. The percentage of cases with Cr* > 1 was 23.41% (1923 out of 8213), which is reasonly in line 
with biomedical observations ||5j[T0l. Similarly, the cases where Cr* < were considered to the ones for 
which sustained virological response cannot be achieved. Samples from each of these sets were compared 
against the frequencies of the biomedically accepted samples. 

We considered values of ep and p in the range [0, 1] and in increments of 0.01. For these pair of 101 x 
101 such values we kept a count of the number of values of sample points which satisfied the condition 
(1 — ep)(l — p) < Cr*. Recall that, this condition corresponds to the stability of the uninfected steady state. 
Thus, this gives us the cumulative distribution of the percentage of cases which give a stable uninfected 
steady state for the corresponding drug efficacies. The results for these are presented in a contour plot in 
Figure ([B and in a surface plot in Figure I©. We plotted the results for the various values of against p 
(Figure O and for values of p against (Figure |4l). One observes that for smaller values of IFN efficacy e^, 
the response (in terms of percentage of cases with stable uninfected steady state) is very sensitive to change 
in the efficacy p of ribavirin. The impact of p however is much less evident in case of IFN being highly 
effective. This observation is consistent with the findings of Dixit et al. El. 



Parameter 


Value na 


s 


1.0 cell mr^ day"\ 


d 


0.01 day-\ 


V 


2.9 virions day~^. 


P 


2.25 X 10"'' ml day~^ virions"^. 


c 


6.0 day-^ 


6 


1.0 day-^ 


r 


2.0 day-i 


T 

^ max 


3.6 X 10^ cells ml-\ 



Table 1 : Mean parameter values for the Latin hypercube sampling 



For the samples which have Cr* values greater than 1 or between and 1, we observe that for variables 
p, (3, c, 6 and Tyaa.x, the null hypothesis (independence) is rejected. It implies that the sample points from 
these two sets tend to deviate more from their original distribution which leads them to be more restricted 
in nature. Sample points which have a negative Cr* value predominantly depend upon the death rates d 
and 6 of the uninfected and infected hepatocytes respectively, for which the null hypothesis fails. Clearly, 
these are cases exhibiting larger values of d and smaller values of 6. In biological terms, these cases present 
with a higher clearance rate of uninfected hepatocytes and a smaller rate of infected hepatocytes. It can 
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also be inferred from equation (12.11 ) that there is little role for drug efficacies if the above conditions are 
prevalent. Also, 6 does not follow the original distribution in any of the cases, indicating its importance in 
the distinction amongst the three cases. 

Next, we chose a threshold value < Crt^r < 1 and analyzed the independence of different parameters 
which have Cr* values on both sides of Crthr- Notice that if Cr* > Crthr, the uninfected steady state can 
be reached for smaller values of drug efficacies as compared to Cr* < Crthr- By varying the value of Crthr 
from 0.2 to 0.4, a Spearman's test showed a modest monotonic relationship between p and /3 for the samples 
with Cr* > Crthr and not in case of Cr* < Crthr- Biologically, /3 represents the rate of infection and p 
represents the rate with which infected hepatocytes are converted into virions. The results highlight the fact 
that for the uninfected steady state to be easily reachable, these two rates cannot be simultaneously high. 

5 Conclusion 

In Latin hypercube sampling, one is at the liberty to choose distributions which best fits the needs (such 
as previous findings, physiological constraints etc.). In this paper, we applied the Latin hypercube sampling 
to generate a large number of sample parameter values for a HCV model. Once the sampling was done, we 
chose values that were positive and satisfied the physiological conditions. One may perform goodness of fit 
and independence tests at every stage of adding constraints to test how much and in what way they affect the 
distribution of the generated samples. This lets us better understand the implications of said constraints on 
a large population of sample parameters. In our case, the values found feasible were subject to the x^-test 
and the Spearman's test. The tests, were very useful in checking the goodness of fit and independence 
of a parameter against an intended distribution. On the other hand. Spearman's test was used to check for 
a monotonic relationship between parameters. The methods and tests mentioned in this paper which were 
also performed on a sample generated from the log-normal distribution gave similar results. 
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Contour plot of percentage of cases with stable uninfected steady state 




Figure 1 : Contour plot of percentage of cases with stable uninfected steady state 
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Percentage of cases with stable uninfected steady state for various e 




Figure 3: Percentage of cases with stable uninfected steady state for various e 
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